This is the analysis for an online survey to gather information about the acceptability of form factor and placements of light loggers. The survey was conducted using Prolific on three days in 2023, each collecting from a different sample (USA, UK, Worldwide).
This is the second part of the analysis, focusing on inference.
We will start off by importing the data and setting factor levels.
#Load Hmisc library
library(Hmisc)
library(tidyverse)
library(gtsummary)
library(gt)
library(ggtext)
library(rlang)
library(patchwork)
library(ordinal)
library(rlang)
library(glue)
library(broom)
library(ggridges)
#Read Data in the long format
load(file = "data/cleaned/LightLoggerFormFacto_Long_CLEANED_DATA_2024-07-02.RData")The base hypothesis is that the ratings depend on the Wearing Position, with Sex and Sample as control variables. Further, the effect of the Wearing Position might be dependent on Sex.
The inferential analysis fits cumulative link mixed models (CLMM) to the data. The models are fitted using the ordinal package. The models are fitted for each rating separately. The base model is fitted with the following formula:
CLMM: \(E(logit(y)) = \alpha_0 + \beta_{i1} * WearingPosition + \beta_{i2} * Sex + \beta_{i3} * WearingPosition*Sex + \beta_{i4} * Sample + \alpha_j\)
Where \(E(logit(y))\) is the expected value of the rating on a logit scale, \(\alpha_0\) is the intercept, \(\beta_{i1}\) is the effect size of Wearing Position, \(\beta_{i2}\) is the effect size of Sex, \(\beta_{i3}\) is the interaction effect size of Wearing Position and Sex, \(\beta_{i4}\) is the effect size of Sample, and \(\alpha_j\) is the random intercept for each participant.
or in Wilkinson notation:
\(Rating \sim WearingPosition*Sex + Sample + (1|Id)\) (m1)
Models are compared through likelihood ratio tests. The following models are used for comparison:
A model without the interaction of Wearing Position and Sex (m2), a model without Sex (m3), a model without Sample (m4), a model with only Wearing Position (m5), and a null model, with only random effects (m0).
\(Rating \sim WearingPosition + Sex + Sample + (1|Id)\) (m2)
\(Rating \sim WearingPosition + Sample + (1|Id)\) (m3)
\(Rating \sim WearingPosition*Sex + (1|Id)\) (m4)
\(Rating \sim WearingPosition + (1|Id)\) (m5)
\(Rating \sim 1 + (1|Id)\) (m0)
The following comparisons are made: 1. Interaction of Sex and Wearing Position (m1 vs. m2) 2. Inclusion of Sex (m2 vs. m3) 3. Inclusion of Sample (m1 vs. m4) 4. Inclusion of Wearing Position (m5 vs. m0)
The p-values are adjusted using the Benjamini & Hochberg (1995) method, also called the false discovery rate (fdr) with an n=4 for four comparisons.
#convert some variables to factor
data_long$sample_location <- as.factor(data_long$sample_location)
data_long$record_id <- as.factor(data_long$record_id)
#function generation
Infe <- function(parameter, data) {
parameter <- enexpr(parameter)
#formula generation
formulas <- inject(c(
#full formula
m1 = !!parameter ~ wearing_position*demographics_sex + sample_location + (1|record_id),
#excluding interaction from m1
m2 = !!parameter ~ wearing_position + demographics_sex + sample_location + (1|record_id),
#further excluding sex from m2
m3 = !!parameter ~ wearing_position + sample_location + (1|record_id),
#excluding sample location from m1
m4 = !!parameter ~ wearing_position*demographics_sex + (1|record_id),
#excluding sample location from m3
m5 = !!parameter ~ wearing_position + (1|record_id),
#Null model
m0 = !!parameter ~ 1 + (1|record_id)
))
#Model generation
Infe_table <- tibble(
name = names(formulas),
formula = formulas,
model = map(name, \(x) {switch(x,
m1 = ,
m2 = ,
m3 = ,
m4 =
clmm(formulas[[x]], data = data, nAGQ = 10,
subset =
demographics_sex != "Other"),
m5 = ,
m0 = clmm(formulas[[x]], data = data, nAGQ = 10)
)
}
)
)
Infe_table
}
# Function to perform anova on model pairs
anova_models <- function(pair) {
one <- pair[[1]]
two <- pair[[2]]
inject(anova(Models$model[[!!one]], Models$model[[!!two]]))
}
#function to perform the model comparisons
Model_Comparisons <- function(Models) {
# Define the model comparisons
model_comparisons <- list(
c(1, 2),
c(2, 3),
c(1, 4),
c(5, 6)
)
tibble(
desc = c("interaction of sex and wearing position",
"inclusion of sex",
"inclusion of sample location",
"inclusion of wearing position"),
p.value = map(model_comparisons,
\(x) anova_models(x) %>% .$`Pr(>Chisq)` %>% .[[2]]
) %>% unlist() %>% p.adjust(method = "fdr", n = 4) %>%
{case_when(
. < 0.001 ~ "<0.001",
TRUE ~ as.character(round(., 3))
)}
)
}
#probability of different answers
#setting up a function
pred <- function(eta, theta, levels, cat = 1:(length(theta)+1), inv.link = plogis) {
Theta <- c(-Inf, theta, Inf)
eta <- c(reference = 0, eta)
table <-
sapply(cat, function(j)
inv.link(Theta[j+1] - eta) - inv.link(Theta[j] - eta) )
colnames(table) <- levels
table
}
data_significance_matrix <- function(Model, subset = "none") {
Models <- tibble(
reference = levels(data_long$wearing_position),
data = map(reference,
\(x) data_long %>% mutate(wearing_position =
fct_relevel(wearing_position, x))
),
model = map(data,
\(x) switch(subset,
none = clmm(Model$formula, data = x, nAGQ = 10),
clmm(Model$formula, data = x, nAGQ = 10,
subset = demographics_sex != "Other")
)
),
tidy = map(model,
\(x) tidy(x) %>%
filter(str_detect(term, "wearing_position")) %>%
select(term, p.value) %>%
mutate(term = str_remove(term, "wearing_position"))
)
)
Models %>%
select(-data, -model) %>%
unnest(tidy) %>%
rowwise() %>%
mutate(
p.value = p.adjust(p.value, method = "fdr", n = 7),
different = p.value <= 0.05) %>%
ungroup() %>%
rbind(
tibble(reference = levels(data_long$wearing_position),
term = reference,
p.value = NA,
different = FALSE))
}
#McFadden´s Pseudo R2, values between 0.2 bis 0.4 are an excellent model fit
McFadden <- function(model, nullmodel) {
ll.null <- nullmodel$logLik
ll.proposed <- model$logLik
1 - ll.proposed/ll.null
}
#Formula to plot the significance matrix
Significance_matrix <- function(sig_matrix, sex = FALSE){
sig_matrix %>%
mutate(
across(c(reference, term), \(x) fct_relevel(x, levels(data_long$wearing_position)))) %>%
ggplot(aes(x=reference, y = term, fill = different)) +
geom_tile(col = "grey30") +
scale_fill_manual(values = c("white", "skyblue4"))+
coord_fixed() +
theme_void() +
labs(title = label(data_long[parameter]),
subtitle = case_when(!sex ~ "<span style='color:skyblue4'>Significant differences</span> due to wearing position",
sex ~ "Significant differences due to <span style='color:skyblue4'>wearing position</span> and <span style='color:red'>sex</span>" )
)+
theme(legend.position = "none",
axis.text.x = element_text(angle = 90, hjust = 1),
axis.text.y = element_text(hjust = 1),
plot.title =
element_textbox_simple(
hjust = 0,
padding = margin(0.5, 0, 0.5, 0, "cm")),
plot.subtitle = element_textbox_simple(),
plot.title.position = "plot",
)
}
#table of probabilities
prob_table <- function(beta, title){
pred(beta, Model$Theta,
levels = levels(data_long[[parameter]])) %>%
as_tibble(rownames = "Wearing Position") %>%
mutate(`Wearing Position` =
str_remove(`Wearing Position`, "wearing_position"),
`Wearing Position` = str_replace(`Wearing Position`, "reference", "Chest Pin")
) %>%
gt() %>% fmt_percent() %>%
gt::tab_caption(md(glue(title)))
}parameter <- "appearance"
parameter_symbol <- sym(parameter)
#model generation
Models <- Infe(!!parameter_symbol, data_long)
#model comparisons
Model_Comparisons(Models) %>%
gt(caption = md(glue("Model comparisons for **{parameter}**"))) %>%
fmt_number(decimals = 3)| desc | p.value |
|---|---|
| interaction of sex and wearing position | 0.827 |
| inclusion of sex | 0.827 |
| inclusion of sample location | 0.57 |
| inclusion of wearing position | <0.001 |
The only parameter which is supported by the data is Wearing Position. Sex, its interaction with Wearing Position, and Sample are not supported by the data. The resulting model m5 is used for the following analysis.
#final model
chosen_model <- "m5"
Model <- Models$model[Models$name == chosen_model][[1]]
summary(Model)Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite
quadrature approximation with 10 quadrature points
formula: appearance ~ wearing_position + (1 | record_id)
data: data
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 1160 -1774.52 3577.05 1604(7930) 2.38e-03 1.6e+02
Random effects:
Groups Name Variance Std.Dev.
record_id (Intercept) 1.531 1.237
Number of groups: record_id 145
Coefficients:
Estimate Std. Error z value Pr(>|z|)
wearing_positionWrist -0.4594 0.2146 -2.141 0.0323 *
wearing_positionNecklace -1.0217 0.2161 -4.728 2.26e-06 ***
wearing_positionSleeve collar -1.5935 0.2253 -7.072 1.53e-12 ***
wearing_positionCollar pin -2.2585 0.2212 -10.211 < 2e-16 ***
wearing_positionNeck loop -2.2304 0.2249 -9.919 < 2e-16 ***
wearing_positionHat pin -2.5887 0.2239 -11.563 < 2e-16 ***
wearing_positionGlasses -3.1912 0.2349 -13.584 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Threshold coefficients:
Estimate Std. Error z value
Extremely unattractive|Very unattractive -4.6499 0.2330 -19.957
Very unattractive|Unattractive -3.7687 0.2179 -17.295
Unattractive|Neutral -2.3197 0.2008 -11.553
Neutral|Attractive -0.6844 0.1886 -3.629
Attractive|Very attractive 1.7319 0.2029 8.537
Very attractive|Extremely attractive 3.3479 0.2744 12.202#McFadden
McFadden(Model, Models$model[[6]])[1] 0.08141085#Significance matrix
sig_matrix <- data_significance_matrix(Model)
P1 <- Significance_matrix(sig_matrix)
ggsave(glue("output/04_unused/significance_matrix_{parameter}.pdf"),
P1, width = 4, height = 4.6)
P1
#if one wanted to print predictions for different subject percentiles
stDev <- as.numeric(attr(VarCorr(Model)$record_id, "stddev"))^2
# 5th percentile:
prob_table(Model$beta * qnorm(0.05) * stDev,
"Predicted probabilities for the **5th percentile** for **{parameter}**")| Wearing Position | Extremely unattractive | Very unattractive | Unattractive | Neutral | Attractive | Very attractive | Extremely attractive |
|---|---|---|---|---|---|---|---|
| Chest Pin | 0.95% | 1.31% | 6.69% | 24.58% | 51.44% | 11.64% | 3.40% |
| Wrist | 0.30% | 0.42% | 2.28% | 10.69% | 50.30% | 25.95% | 10.05% |
| Necklace | 0.07% | 0.10% | 0.57% | 2.96% | 26.44% | 38.33% | 31.52% |
| Sleeve collar | 0.02% | 0.02% | 0.14% | 0.73% | 8.38% | 24.71% | 66.01% |
| Collar pin | 0.00% | 0.00% | 0.03% | 0.14% | 1.71% | 6.92% | 91.20% |
| Neck loop | 0.00% | 0.00% | 0.03% | 0.15% | 1.83% | 7.37% | 90.61% |
| Hat pin | 0.00% | 0.00% | 0.01% | 0.06% | 0.75% | 3.21% | 95.97% |
| Glasses | 0.00% | 0.00% | 0.00% | 0.01% | 0.17% | 0.73% | 99.09% |
#average:
prob_table(Model$beta,
"Predicted probabilities for the **average person** for **{parameter}**")| Wearing Position | Extremely unattractive | Very unattractive | Unattractive | Neutral | Attractive | Very attractive | Extremely attractive |
|---|---|---|---|---|---|---|---|
| Chest Pin | 0.95% | 1.31% | 6.69% | 24.58% | 51.44% | 11.64% | 3.40% |
| Wrist | 1.49% | 2.03% | 9.94% | 30.93% | 45.55% | 7.88% | 2.17% |
| Necklace | 2.59% | 3.44% | 15.42% | 36.90% | 35.66% | 4.74% | 1.25% |
| Sleeve collar | 4.49% | 5.71% | 22.40% | 38.68% | 25.25% | 2.76% | 0.71% |
| Collar pin | 8.38% | 9.71% | 30.38% | 34.37% | 15.35% | 1.45% | 0.37% |
| Neck loop | 8.17% | 9.51% | 30.09% | 34.66% | 15.70% | 1.49% | 0.38% |
| Hat pin | 11.29% | 12.21% | 33.18% | 30.35% | 11.65% | 1.05% | 0.26% |
| Glasses | 18.87% | 17.08% | 34.55% | 21.96% | 6.82% | 0.58% | 0.14% |
# 95th percentile:
prob_table(Model$beta * qnorm(0.95) * stDev,
"Predicted probabilities for the **95th percentile** for **{parameter}**")| Wearing Position | Extremely unattractive | Very unattractive | Unattractive | Neutral | Attractive | Very attractive | Extremely attractive |
|---|---|---|---|---|---|---|---|
| Chest Pin | 0.95% | 1.31% | 6.69% | 24.58% | 51.44% | 11.64% | 3.40% |
| Wrist | 2.95% | 3.89% | 16.97% | 37.78% | 33.14% | 4.18% | 1.09% |
| Necklace | 11.13% | 12.08% | 33.07% | 30.57% | 11.82% | 1.07% | 0.27% |
| Sleeve collar | 34.56% | 21.48% | 28.40% | 12.09% | 3.15% | 0.26% | 0.06% |
| Collar pin | 73.80% | 13.38% | 9.48% | 2.67% | 0.61% | 0.05% | 0.01% |
| Neck loop | 72.41% | 13.96% | 10.06% | 2.86% | 0.65% | 0.05% | 0.01% |
| Hat pin | 86.61% | 7.37% | 4.54% | 1.19% | 0.27% | 0.02% | 0.01% |
| Glasses | 96.72% | 1.89% | 1.06% | 0.26% | 0.06% | 0.00% | 0.00% |
#random effect distribution with random effects
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:145, labels=ord.re)
axis(2)
for(i in 1:145) segments(i, ci[i,1], i, ci[i,2])
abline(h=0, lty=2)
#without confidence intervals
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:145, labels=ord.re)
axis(2)
abline(h=0, lty=2)
This sections analyses how the context of wear influences the ratings, dependent on the wearing position
#function to remove labels. This is necessary for pivot_longer to work
clear.labels <- function(x) {
if(is.list(x)) {
for(i in seq_along(x)) {
class(x[[i]]) <- setdiff(class(x[[i]]), 'labelled')
attr(x[[i]],"label") <- NULL
}
} else {
class(x) <- setdiff(class(x), "labelled")
attr(x, "label") <- NULL
}
return(x)
}
#pivoting the data even longer to get the context into the observation
data_long2 <- data_long %>%
select(record_id,
demographics_sex,
wearing_position,
starts_with("wearlikelihoodrating")
) %>%
clear.labels() %>%
pivot_longer(-c(record_id, demographics_sex, wearing_position),
names_to = "context",
values_to = "rating") %>%
mutate(context =
str_remove(context, "wearlikelihoodrating_") %>%
factor() %>%
fct_reorder(as.numeric(rating), .desc = TRUE, .fun = mean)
# fct_relevel(c("public", "work", "home", "social", "exercise")
)
#model generation
Model <- clmm(rating ~ context*wearing_position + (1|record_id), data = data_long2, nAGQ = 10)
Model2 <- clmm(rating ~ context+wearing_position + (1|record_id), data = data_long2, nAGQ = 10)
Model3 <- clmm(rating ~ context + (1|record_id), data = data_long2, nAGQ = 10)
Model4 <- clmm(rating ~ wearing_position + (1|record_id), data = data_long2, nAGQ = 10)
Model0 <- clmm(rating ~ 1 + (1|record_id), data = data_long2, nAGQ = 10)
#model selection
anova(Model, Model2)Likelihood ratio tests of cumulative link models:
formula: link: threshold:
Model2 rating ~ context + wearing_position + (1 | record_id) logit flexible
Model rating ~ context * wearing_position + (1 | record_id) logit flexible
no.par AIC logLik LR.stat df Pr(>Chisq)
Model2 18 18590 -9276.8
Model 46 18448 -9178.1 197.41 28 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(Model2, Model3)Likelihood ratio tests of cumulative link models:
formula: link: threshold:
Model3 rating ~ context + (1 | record_id) logit flexible
Model2 rating ~ context + wearing_position + (1 | record_id) logit flexible
no.par AIC logLik LR.stat df Pr(>Chisq)
Model3 11 19937 -9957.5
Model2 18 18590 -9276.8 1361.5 7 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(Model3, Model0)Likelihood ratio tests of cumulative link models:
formula: link: threshold:
Model0 rating ~ 1 + (1 | record_id) logit flexible
Model3 rating ~ context + (1 | record_id) logit flexible
no.par AIC logLik LR.stat df Pr(>Chisq)
Model0 7 20120 -10053.2
Model3 11 19937 -9957.5 191.28 4 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(Model4, Model0)Likelihood ratio tests of cumulative link models:
formula: link: threshold:
Model0 rating ~ 1 + (1 | record_id) logit flexible
Model4 rating ~ wearing_position + (1 | record_id) logit flexible
no.par AIC logLik LR.stat df Pr(>Chisq)
Model0 7 20120 -10053.2
Model4 14 18831 -9401.6 1303.1 7 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#model summaries
summary(Model)Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite
quadrature approximation with 10 quadrature points
formula: rating ~ context * wearing_position + (1 | record_id)
data: data_long2
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 5800 -9178.09 18448.18 12392(93749) 6.64e-03 3.9e+03
Random effects:
Groups Name Variance Std.Dev.
record_id (Intercept) 2.026 1.423
Number of groups: record_id 145
Coefficients:
Estimate Std. Error z value
contextpublic -0.4639 0.2067 -2.244
contextexercise -0.5336 0.2125 -2.511
contextwork -0.7212 0.2086 -3.458
contextsocial -0.9508 0.2078 -4.576
wearing_positionWrist -0.7158 0.2148 -3.333
wearing_positionNecklace -0.6401 0.2137 -2.996
wearing_positionSleeve collar -1.4403 0.2170 -6.636
wearing_positionCollar pin -1.3705 0.2114 -6.481
wearing_positionNeck loop -1.7915 0.2148 -8.338
wearing_positionHat pin -2.9718 0.2167 -13.713
wearing_positionGlasses -2.4774 0.2142 -11.569
contextpublic:wearing_positionWrist 0.3498 0.2977 1.175
contextexercise:wearing_positionWrist 0.8694 0.3055 2.846
contextwork:wearing_positionWrist 0.2889 0.3021 0.956
contextsocial:wearing_positionWrist 0.3082 0.2979 1.035
contextpublic:wearing_positionNecklace -0.1876 0.2965 -0.633
contextexercise:wearing_positionNecklace -1.3970 0.3048 -4.584
contextwork:wearing_positionNecklace -0.3569 0.2996 -1.191
contextsocial:wearing_positionNecklace -0.2508 0.2978 -0.842
contextpublic:wearing_positionSleeve collar -0.2226 0.3019 -0.737
contextexercise:wearing_positionSleeve collar 1.3441 0.3107 4.326
contextwork:wearing_positionSleeve collar -0.5646 0.3007 -1.877
contextsocial:wearing_positionSleeve collar -0.4524 0.3004 -1.506
contextpublic:wearing_positionCollar pin -0.5273 0.2931 -1.799
contextexercise:wearing_positionCollar pin 0.1147 0.3011 0.381
contextwork:wearing_positionCollar pin -0.6189 0.2981 -2.076
contextsocial:wearing_positionCollar pin -0.6357 0.2964 -2.145
contextpublic:wearing_positionNeck loop -0.2973 0.2964 -1.003
contextexercise:wearing_positionNeck loop -0.3186 0.3064 -1.040
contextwork:wearing_positionNeck loop -0.5299 0.3005 -1.763
contextsocial:wearing_positionNeck loop -0.3003 0.2990 -1.005
contextpublic:wearing_positionHat pin 0.5526 0.2968 1.862
contextexercise:wearing_positionHat pin 0.3640 0.3055 1.192
contextwork:wearing_positionHat pin -0.1569 0.3023 -0.519
contextsocial:wearing_positionHat pin 0.4544 0.3003 1.513
contextpublic:wearing_positionGlasses -0.2479 0.2995 -0.828
contextexercise:wearing_positionGlasses -0.7585 0.3049 -2.488
contextwork:wearing_positionGlasses -0.1357 0.3020 -0.449
contextsocial:wearing_positionGlasses -0.2515 0.3036 -0.829
Pr(>|z|)
contextpublic 0.024811 *
contextexercise 0.012028 *
contextwork 0.000545 ***
contextsocial 4.74e-06 ***
wearing_positionWrist 0.000860 ***
wearing_positionNecklace 0.002735 **
wearing_positionSleeve collar 3.21e-11 ***
wearing_positionCollar pin 9.09e-11 ***
wearing_positionNeck loop < 2e-16 ***
wearing_positionHat pin < 2e-16 ***
wearing_positionGlasses < 2e-16 ***
contextpublic:wearing_positionWrist 0.239885
contextexercise:wearing_positionWrist 0.004433 **
contextwork:wearing_positionWrist 0.338828
contextsocial:wearing_positionWrist 0.300884
contextpublic:wearing_positionNecklace 0.526902
contextexercise:wearing_positionNecklace 4.56e-06 ***
contextwork:wearing_positionNecklace 0.233557
contextsocial:wearing_positionNecklace 0.399794
contextpublic:wearing_positionSleeve collar 0.460819
contextexercise:wearing_positionSleeve collar 1.52e-05 ***
contextwork:wearing_positionSleeve collar 0.060497 .
contextsocial:wearing_positionSleeve collar 0.132102
contextpublic:wearing_positionCollar pin 0.072055 .
contextexercise:wearing_positionCollar pin 0.703195
contextwork:wearing_positionCollar pin 0.037868 *
contextsocial:wearing_positionCollar pin 0.031984 *
contextpublic:wearing_positionNeck loop 0.315737
contextexercise:wearing_positionNeck loop 0.298500
contextwork:wearing_positionNeck loop 0.077833 .
contextsocial:wearing_positionNeck loop 0.315105
contextpublic:wearing_positionHat pin 0.062620 .
contextexercise:wearing_positionHat pin 0.233364
contextwork:wearing_positionHat pin 0.603725
contextsocial:wearing_positionHat pin 0.130261
contextpublic:wearing_positionGlasses 0.407884
contextexercise:wearing_positionGlasses 0.012838 *
contextwork:wearing_positionGlasses 0.653081
contextsocial:wearing_positionGlasses 0.407351
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Threshold coefficients:
Estimate Std. Error z value
Extremely unlikely|Very Unlikely -4.2432 0.1984 -21.386
Very Unlikely|Unlikely -3.3268 0.1960 -16.973
Unlikely|Neutral -2.2250 0.1940 -11.470
Neutral|Likely -1.2738 0.1927 -6.609
Likely|Very likely 0.4591 0.1917 2.395
Very likely|Extremely likely 2.1776 0.1976 11.021summary(Model3)Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite
quadrature approximation with 10 quadrature points
formula: rating ~ context + (1 | record_id)
data: data_long2
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 5800 -9957.53 19937.05 1209(8342) 8.62e-03 3.6e+02
Random effects:
Groups Name Variance Std.Dev.
record_id (Intercept) 1.451 1.205
Number of groups: record_id 145
Coefficients:
Estimate Std. Error z value Pr(>|z|)
contextpublic -0.45146 0.07412 -6.091 1.12e-09 ***
contextexercise -0.44588 0.07606 -5.862 4.57e-09 ***
contextwork -0.83751 0.07543 -11.103 < 2e-16 ***
contextsocial -0.91921 0.07517 -12.229 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Threshold coefficients:
Estimate Std. Error z value
Extremely unlikely|Very Unlikely -2.3870 0.1188 -20.100
Very Unlikely|Unlikely -1.5906 0.1163 -13.674
Unlikely|Neutral -0.6642 0.1149 -5.782
Neutral|Likely 0.1109 0.1145 0.968
Likely|Very likely 1.5633 0.1167 13.391
Very likely|Extremely likely 3.1163 0.1288 24.187summary(Model4)Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite
quadrature approximation with 10 quadrature points
formula: rating ~ wearing_position + (1 | record_id)
data: data_long2
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 5800 -9401.64 18831.28 2076(16279) 9.58e-03 4.0e+02
Random effects:
Groups Name Variance Std.Dev.
record_id (Intercept) 1.851 1.36
Number of groups: record_id 145
Coefficients:
Estimate Std. Error z value Pr(>|z|)
wearing_positionWrist -0.33834 0.09403 -3.598 0.000321 ***
wearing_positionNecklace -1.01892 0.09474 -10.755 < 2e-16 ***
wearing_positionSleeve collar -1.39982 0.09718 -14.405 < 2e-16 ***
wearing_positionCollar pin -1.62095 0.09572 -16.935 < 2e-16 ***
wearing_positionNeck loop -1.99074 0.09701 -20.522 < 2e-16 ***
wearing_positionHat pin -2.59863 0.09921 -26.194 < 2e-16 ***
wearing_positionGlasses -2.62203 0.10032 -26.136 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Threshold coefficients:
Estimate Std. Error z value
Extremely unlikely|Very Unlikely -3.5293 0.1392 -25.352
Very Unlikely|Unlikely -2.6479 0.1362 -19.441
Unlikely|Neutral -1.6022 0.1338 -11.972
Neutral|Likely -0.7050 0.1324 -5.323
Likely|Very likely 0.9436 0.1327 7.109
Very likely|Extremely likely 2.6047 0.1431 18.196#get the predicted probabilities depending on context, save as a table
context_preds <- pred(Model3$beta, Model3$Theta,
levels = levels(data_long2$rating)) %>%
as_tibble(rownames = "Parameters") %>%
mutate(`Parameters` =
str_remove(`Parameters`, "wearing_position|context"),
`Parameters` = str_replace(`Parameters`, "reference", "home")
)
context_preds_table <-
context_preds %>%
gt() %>%
fmt_percent(decimals = 0) %>%
gt::tab_caption(md("Likelihood of wear ratings depending on context")) %>%
gt::data_color(columns = -Parameters, method = "numeric",
colors = scales::col_numeric(palette = "viridis",
domain = c(0,0.30))
) %>%
cols_width(everything() ~ px(100))Warning: Since gt v0.9.0, the `colors` argument has been deprecated.
• Please use the `fn` argument instead.
This warning is displayed once every 8 hours.context_preds_table| Parameters | Extremely unlikely | Very Unlikely | Unlikely | Neutral | Likely | Very likely | Extremely likely |
|---|---|---|---|---|---|---|---|
| home | 8% | 9% | 17% | 19% | 30% | 13% | 4% |
| public | 13% | 12% | 20% | 19% | 25% | 9% | 3% |
| exercise | 13% | 12% | 20% | 19% | 25% | 9% | 3% |
| work | 18% | 14% | 22% | 18% | 20% | 6% | 2% |
| social | 19% | 15% | 23% | 17% | 19% | 6% | 2% |
gtsave(context_preds_table, "output/02_tables/table_2_context.png")file:////var/folders/9p/326_k3kx43qbn_cyl1rqfhb00000gn/T//Rtmp2v7SBD/file85a86024d65a.html screenshot completed#get the predicted probabilities depending on wearing_position, save as a table
context_preds <- pred(Model4$beta, Model4$Theta,
levels = levels(data_long2$rating)) %>%
as_tibble(rownames = "Parameters") %>%
mutate(`Parameters` =
str_remove(`Parameters`, "wearing_position|context"),
`Parameters` = str_replace(`Parameters`, "reference", "Chest pin")
)
context_preds_table <-
context_preds %>%
gt() %>%
fmt_percent(decimals = 0) %>%
gt::tab_caption(md("Likelihood of wear ratings depending on wearing position")) %>%
gt::data_color(columns = -Parameters, method = "numeric",
colors = scales::col_numeric(palette = "viridis",
domain = c(0,0.39))
) %>%
cols_width(everything() ~ px(100))
context_preds_table| Parameters | Extremely unlikely | Very Unlikely | Unlikely | Neutral | Likely | Very likely | Extremely likely |
|---|---|---|---|---|---|---|---|
| Chest pin | 3% | 4% | 10% | 16% | 39% | 21% | 7% |
| Wrist | 4% | 5% | 13% | 19% | 37% | 17% | 5% |
| Necklace | 8% | 9% | 19% | 22% | 30% | 10% | 3% |
| Sleeve collar | 11% | 12% | 23% | 22% | 25% | 7% | 2% |
| Collar pin | 13% | 13% | 24% | 21% | 21% | 6% | 1% |
| Neck loop | 18% | 16% | 25% | 19% | 17% | 4% | 1% |
| Hat pin | 28% | 20% | 24% | 14% | 10% | 2% | 1% |
| Glasses | 29% | 21% | 24% | 14% | 10% | 2% | 1% |
gtsave(context_preds_table, "output/02_tables/table_3_wearing_position.png")file:////var/folders/9p/326_k3kx43qbn_cyl1rqfhb00000gn/T//Rtmp2v7SBD/file85a86e794f42.html screenshot completed#create a table of the interaction model
#first, collect all the beta coefficients for the context
coefs_context <-
Model$beta %>%
as_tibble(rownames = "Parameters") %>%
filter(str_detect(Parameters, "context"),
!str_detect(Parameters, "wearing_position")) %>%
rename(Context = Parameters,
beta1 = value) %>%
add_row(Context = "contexthome", beta1 = 0, .before = 1)
#second, collect all the beta coefficients for the wearing position
coefs_position <-
Model$beta %>%
as_tibble(rownames = "Parameters") %>%
filter(str_detect(Parameters, "wearing_position"),
!str_detect(Parameters, "context")) %>%
rename(Position2 = Parameters,
beta2 = value) %>%
add_row(Position2 = "wearing_positionChest pin", beta2 = 0, .before = 1)
#third, collect all the beta coefficients for the interaction
coefs_extra <-
Model$beta %>%
as_tibble(rownames = "Parameters") %>%
filter(str_detect(Parameters, "wearing_position")) %>%
rename(Position = Parameters,
beta3 = value) %>%
add_row(Position = "contextpublic:wearing_positionChest pin", beta3 = 0, .before = 2) %>%
add_row(Position = "contextexercise:wearing_positionChest pin", beta3 = 0, .before = 3) %>%
add_row(Position = "contextwork:wearing_positionChest pin", beta3 = 0, .before = 4) %>%
add_row(Position = "contextsocial:wearing_positionChest pin", beta3 = 0, .before = 5) %>%
mutate(Position =
str_replace(
Position, "^wearing_position", "contexthome:wearing_position")
)
#combine the three tables to a single table and filter out non-relevant combinations
coefs_combined <-
expand_grid(coefs_context, coefs_position, coefs_extra) %>%
filter(str_detect(Position, Context),
str_detect(Position, Position2)) %>%
mutate(Context = str_remove(Context, "context"),
Position = str_remove(Position, "(.*):wearing_position"),
beta3 = ifelse(str_detect(Context , "home"), 0, beta3))
coefs_combined <-
set_names(coefs_combined$beta1+coefs_combined$beta2 + coefs_combined$beta3,
paste(coefs_combined$Context, coefs_combined$Position, sep = ":"))
#create and save a table of the predicted rating
combined_pred <-
pred(coefs_combined, Model$Theta,
levels = levels(data_long2$rating)
) %>%
as_tibble(rownames = "Parameters") %>%
mutate(
`Parameters` = str_replace(`Parameters`, "reference", "home:Chest pin")
) %>%
separate_wider_delim(Parameters, ":", names = c("Context", "Position")) %>%
pivot_longer(cols = -c(Context, Position), names_to = "Rating", values_to = "Probability") %>%
mutate(Rating = factor(Rating, levels = levels(data_long2$rating))) %>%
group_by(Context, Position) %>%
summarize(Probability = Rating[which.max(Probability)], .groups = "drop") %>%
pivot_wider(names_from = Position, values_from = Probability) %>%
arrange(factor(Context, levels = levels(data_long2$context))) %>%
select(Context, all_of(levels(data_long2$wearing_position))) %>%
gt() %>%
gt::tab_caption(md("*Likelihood of wear* rating prediction depending on **context** and **wearing position**")) %>%
gt::data_color(columns = -Context, ordered = TRUE,
palette = "viridis"
) %>%
cols_width(everything() ~ px(100))
combined_pred| Context | Chest pin | Wrist | Necklace | Sleeve collar | Collar pin | Neck loop | Hat pin | Glasses |
|---|---|---|---|---|---|---|---|---|
| home | Likely | Likely | Likely | Likely | Likely | Likely | Unlikely | Unlikely |
| public | Likely | Likely | Likely | Unlikely | Unlikely | Unlikely | Unlikely | Extremely unlikely |
| exercise | Likely | Likely | Unlikely | Likely | Likely | Unlikely | Unlikely | Extremely unlikely |
| work | Likely | Likely | Likely | Unlikely | Unlikely | Unlikely | Extremely unlikely | Extremely unlikely |
| social | Likely | Likely | Likely | Unlikely | Unlikely | Unlikely | Extremely unlikely | Extremely unlikely |
gtsave(combined_pred, "output/02_tables/table_4_context_position.png")file:////var/folders/9p/326_k3kx43qbn_cyl1rqfhb00000gn/T//Rtmp2v7SBD/file85a85e42aeb7.html screenshot completedparameter <- "wearlikelihoodrating_public"
parameter_symbol <- sym(parameter)
#model generation
Models <- Infe(!!parameter_symbol, data_long)
#model comparisons
Model_Comparisons(Models) %>%
gt(caption = md(glue("Model comparisons for **{parameter}**"))) %>%
fmt_number(decimals = 3)| desc | p.value |
|---|---|
| interaction of sex and wearing position | 0.812 |
| inclusion of sex | 0.812 |
| inclusion of sample location | 0.389 |
| inclusion of wearing position | <0.001 |
The only parameter which is supported by the data is Wearing Position. Sex, its interaction with Wearing Position, and Sample are not supported by the data. The resulting model m5 is used for the following analysis.
#final model
chosen_model <- "m5"
Model <- Models$model[Models$name == chosen_model][[1]]
summary(Model)Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite
quadrature approximation with 10 quadrature points
formula: wearlikelihoodrating_public ~ wearing_position + (1 | record_id)
data: data
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 1160 -1899.19 3826.39 1570(9364) 1.01e-03 2.2e+02
Random effects:
Groups Name Variance Std.Dev.
record_id (Intercept) 1.921 1.386
Number of groups: record_id 145
Coefficients:
Estimate Std. Error z value Pr(>|z|)
wearing_positionWrist -0.3931 0.2132 -1.844 0.0652 .
wearing_positionNecklace -0.8693 0.2138 -4.067 4.77e-05 ***
wearing_positionSleeve collar -1.7525 0.2218 -7.901 2.76e-15 ***
wearing_positionCollar pin -2.0125 0.2176 -9.250 < 2e-16 ***
wearing_positionNeck loop -2.1992 0.2195 -10.021 < 2e-16 ***
wearing_positionHat pin -2.5519 0.2234 -11.425 < 2e-16 ***
wearing_positionGlasses -2.8638 0.2317 -12.358 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Threshold coefficients:
Estimate Std. Error z value
Extremely unlikely|Very Unlikely -4.0037 0.2262 -17.699
Very Unlikely|Unlikely -3.1317 0.2140 -14.633
Unlikely|Neutral -1.9030 0.2022 -9.410
Neutral|Likely -0.9430 0.1958 -4.816
Likely|Very likely 1.0424 0.1964 5.307
Very likely|Extremely likely 2.8424 0.2398 11.855#McFadden
McFadden(Model, Models$model[[6]])[1] 0.07070644#Significance matrix
sig_matrix <- data_significance_matrix(Model)
P2 <- Significance_matrix(sig_matrix)
ggsave(glue("output/04_unused/significance_matrix_{parameter}.pdf"),
P2, width = 4, height = 4.7)
P2
#if one wanted to print predictions for different subject percentiles
stDev <- as.numeric(attr(VarCorr(Model)$record_id, "stddev"))^2
# 5th percentile:
prob_table(Model$beta * qnorm(0.05) * stDev,
"Predicted probabilities for the **5th percentile** for **{parameter}**")| Wearing Position | Extremely unlikely | Very Unlikely | Unlikely | Neutral | Likely | Very likely | Extremely likely |
|---|---|---|---|---|---|---|---|
| Chest Pin | 1.79% | 2.39% | 8.79% | 15.05% | 45.90% | 20.56% | 5.51% |
| Wrist | 0.52% | 0.72% | 2.88% | 5.98% | 34.91% | 38.18% | 16.80% |
| Necklace | 0.12% | 0.16% | 0.67% | 1.49% | 12.95% | 36.99% | 47.63% |
| Sleeve collar | 0.01% | 0.01% | 0.04% | 0.09% | 0.95% | 5.22% | 93.68% |
| Collar pin | 0.00% | 0.00% | 0.02% | 0.04% | 0.42% | 2.39% | 97.12% |
| Neck loop | 0.00% | 0.00% | 0.01% | 0.02% | 0.23% | 1.35% | 98.38% |
| Hat pin | 0.00% | 0.00% | 0.00% | 0.01% | 0.08% | 0.45% | 99.46% |
| Glasses | 0.00% | 0.00% | 0.00% | 0.00% | 0.03% | 0.17% | 99.80% |
#average:
prob_table(Model$beta,
"Predicted probabilities for the **average person** for **{parameter}**")| Wearing Position | Extremely unlikely | Very Unlikely | Unlikely | Neutral | Likely | Very likely | Extremely likely |
|---|---|---|---|---|---|---|---|
| Chest Pin | 1.79% | 2.39% | 8.79% | 15.05% | 45.90% | 20.56% | 5.51% |
| Wrist | 2.63% | 3.44% | 12.02% | 18.49% | 44.19% | 15.44% | 3.79% |
| Necklace | 4.17% | 5.26% | 16.81% | 21.92% | 38.96% | 10.49% | 2.39% |
| Sleeve collar | 9.52% | 10.59% | 26.13% | 22.96% | 25.04% | 4.76% | 1.00% |
| Collar pin | 12.01% | 12.60% | 28.12% | 21.72% | 21.05% | 3.73% | 0.77% |
| Neck loop | 14.13% | 14.11% | 29.11% | 20.49% | 18.40% | 3.12% | 0.64% |
| Hat pin | 18.97% | 16.93% | 29.78% | 17.65% | 14.00% | 2.22% | 0.45% |
| Glasses | 24.23% | 19.11% | 28.99% | 14.89% | 10.80% | 1.64% | 0.33% |
# 95th percentile:
prob_table(Model$beta * qnorm(0.95) * stDev,
"Predicted probabilities for the **95th percentile** for **{parameter}**")| Wearing Position | Extremely unlikely | Very Unlikely | Unlikely | Neutral | Likely | Very likely | Extremely likely |
|---|---|---|---|---|---|---|---|
| Chest Pin | 1.79% | 2.39% | 8.79% | 15.05% | 45.90% | 20.56% | 5.51% |
| Wrist | 5.94% | 7.19% | 20.93% | 23.37% | 33.33% | 7.59% | 1.66% |
| Necklace | 22.16% | 18.35% | 29.43% | 15.93% | 11.92% | 1.84% | 0.37% |
| Sleeve collar | 82.27% | 9.46% | 5.70% | 1.57% | 0.86% | 0.12% | 0.02% |
| Collar pin | 91.35% | 4.84% | 2.66% | 0.70% | 0.38% | 0.05% | 0.01% |
| Neck loop | 95.01% | 2.84% | 1.51% | 0.39% | 0.21% | 0.03% | 0.01% |
| Hat pin | 98.31% | 0.98% | 0.50% | 0.13% | 0.07% | 0.01% | 0.00% |
| Glasses | 99.36% | 0.37% | 0.19% | 0.05% | 0.03% | 0.00% | 0.00% |
#random effect distribution with random effects
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:145, labels=ord.re)
axis(2)
for(i in 1:145) segments(i, ci[i,1], i, ci[i,2])
abline(h=0, lty=2)
#without confidence intervals
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:145, labels=ord.re)
axis(2)
abline(h=0, lty=2)
parameter <- "wearlikelihoodrating_work"
parameter_symbol <- sym(parameter)
#model generation
Models <- Infe(!!parameter_symbol, data_long)
#model comparisons
Model_Comparisons(Models) %>%
gt(caption = md(glue("Model comparisons for **{parameter}**"))) %>%
fmt_number(decimals = 3)| desc | p.value |
|---|---|
| interaction of sex and wearing position | 0.118 |
| inclusion of sex | 0.648 |
| inclusion of sample location | 0.077 |
| inclusion of wearing position | <0.001 |
The only parameter which is supported by the data is Wearing Position. Sex, its interaction with Wearing Position, and Sample are not supported by the data. The resulting model m5 is used for the following analysis.
#final model
chosen_model <- "m5"
Model <- Models$model[Models$name == chosen_model][[1]]
summary(Model)Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite
quadrature approximation with 10 quadrature points
formula: wearlikelihoodrating_work ~ wearing_position + (1 | record_id)
data: data
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 1160 -1874.88 3777.76 1567(7749) 1.18e-03 2.1e+02
Random effects:
Groups Name Variance Std.Dev.
record_id (Intercept) 2.229 1.493
Number of groups: record_id 145
Coefficients:
Estimate Std. Error z value Pr(>|z|)
wearing_positionWrist -0.4633 0.2148 -2.157 0.031 *
wearing_positionNecklace -1.0287 0.2140 -4.806 1.54e-06 ***
wearing_positionSleeve collar -2.0777 0.2201 -9.441 < 2e-16 ***
wearing_positionCollar pin -2.0459 0.2188 -9.351 < 2e-16 ***
wearing_positionNeck loop -2.4086 0.2227 -10.818 < 2e-16 ***
wearing_positionHat pin -3.2622 0.2331 -13.997 < 2e-16 ***
wearing_positionGlasses -2.7678 0.2303 -12.019 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Threshold coefficients:
Estimate Std. Error z value
Extremely unlikely|Very Unlikely -3.7448 0.2260 -16.570
Very Unlikely|Unlikely -2.6696 0.2135 -12.506
Unlikely|Neutral -1.5506 0.2046 -7.577
Neutral|Likely -0.5348 0.1989 -2.688
Likely|Very likely 1.2863 0.2036 6.317
Very likely|Extremely likely 2.9820 0.2503 11.913#McFadden
McFadden(Model, Models$model[[6]])[1] 0.08277394#Significance matrix
sig_matrix <- data_significance_matrix(Model)
P3 <- Significance_matrix(sig_matrix)
ggsave(glue("output/04_unused/significance_matrix_{parameter}.pdf"),
P3, width = 4, height = 4.7)
P3
#if one wanted to print predictions for different subject percentiles
stDev <- as.numeric(attr(VarCorr(Model)$record_id, "stddev"))^2
# 5th percentile:
prob_table(Model$beta * qnorm(0.05) * stDev,
"Predicted probabilities for the **5th percentile** for **{parameter}**")| Wearing Position | Extremely unlikely | Very Unlikely | Unlikely | Neutral | Likely | Very likely | Extremely likely |
|---|---|---|---|---|---|---|---|
| Chest Pin | 2.31% | 4.17% | 11.02% | 19.44% | 41.41% | 16.82% | 4.82% |
| Wrist | 0.43% | 0.82% | 2.48% | 5.94% | 30.15% | 38.47% | 21.70% |
| Necklace | 0.05% | 0.10% | 0.33% | 0.84% | 6.36% | 23.53% | 68.78% |
| Sleeve collar | 0.00% | 0.00% | 0.01% | 0.02% | 0.15% | 0.78% | 99.04% |
| Collar pin | 0.00% | 0.00% | 0.01% | 0.02% | 0.17% | 0.88% | 98.92% |
| Neck loop | 0.00% | 0.00% | 0.00% | 0.01% | 0.04% | 0.23% | 99.71% |
| Hat pin | 0.00% | 0.00% | 0.00% | 0.00% | 0.00% | 0.01% | 99.99% |
| Glasses | 0.00% | 0.00% | 0.00% | 0.00% | 0.01% | 0.06% | 99.92% |
#average:
prob_table(Model$beta,
"Predicted probabilities for the **average person** for **{parameter}**")| Wearing Position | Extremely unlikely | Very Unlikely | Unlikely | Neutral | Likely | Very likely | Extremely likely |
|---|---|---|---|---|---|---|---|
| Chest Pin | 2.31% | 4.17% | 11.02% | 19.44% | 41.41% | 16.82% | 4.82% |
| Wrist | 3.62% | 6.30% | 15.29% | 23.00% | 36.98% | 11.72% | 3.09% |
| Necklace | 6.20% | 10.03% | 21.01% | 24.86% | 28.91% | 7.21% | 1.78% |
| Sleeve collar | 15.88% | 19.74% | 27.26% | 19.51% | 14.27% | 2.71% | 0.63% |
| Collar pin | 15.46% | 19.43% | 27.24% | 19.79% | 14.63% | 2.80% | 0.65% |
| Neck loop | 20.81% | 22.70% | 26.71% | 16.46% | 10.89% | 1.97% | 0.45% |
| Hat pin | 38.16% | 26.23% | 20.31% | 9.16% | 5.09% | 0.85% | 0.19% |
| Glasses | 27.35% | 25.11% | 24.70% | 13.16% | 7.98% | 1.39% | 0.32% |
# 95th percentile:
prob_table(Model$beta * qnorm(0.95) * stDev,
"Predicted probabilities for the **95th percentile** for **{parameter}**")| Wearing Position | Extremely unlikely | Very Unlikely | Unlikely | Neutral | Likely | Very likely | Extremely likely |
|---|---|---|---|---|---|---|---|
| Chest Pin | 2.31% | 4.17% | 11.02% | 19.44% | 41.41% | 16.82% | 4.82% |
| Wrist | 11.45% | 16.03% | 26.23% | 22.51% | 18.98% | 3.89% | 0.92% |
| Necklace | 50.68% | 24.39% | 15.14% | 6.00% | 3.15% | 0.52% | 0.12% |
| Sleeve collar | 97.96% | 1.33% | 0.47% | 0.15% | 0.07% | 0.01% | 0.00% |
| Collar pin | 97.72% | 1.49% | 0.53% | 0.17% | 0.08% | 0.01% | 0.00% |
| Neck loop | 99.39% | 0.40% | 0.14% | 0.04% | 0.02% | 0.00% | 0.00% |
| Hat pin | 99.97% | 0.02% | 0.01% | 0.00% | 0.00% | 0.00% | 0.00% |
| Glasses | 99.83% | 0.11% | 0.04% | 0.01% | 0.01% | 0.00% | 0.00% |
#random effect distribution with random effects
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:145, labels=ord.re)
axis(2)
for(i in 1:145) segments(i, ci[i,1], i, ci[i,2])
abline(h=0, lty=2)
#without confidence intervals
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:145, labels=ord.re)
axis(2)
abline(h=0, lty=2)
parameter <- "wearlikelihoodrating_home"
parameter_symbol <- sym(parameter)
#model generation
Models <- Infe(!!parameter_symbol, data_long)
#model comparisons
Model_Comparisons(Models) %>%
gt(caption = md(glue("Model comparisons for **{parameter}**"))) %>%
fmt_number(decimals = 3)| desc | p.value |
|---|---|
| interaction of sex and wearing position | 0.332 |
| inclusion of sex | 0.653 |
| inclusion of sample location | 0.943 |
| inclusion of wearing position | <0.001 |
The only parameter which is supported by the data is Wearing Position. Sex, its interaction with Wearing Position, and Sample are not supported by the data. The resulting model m5 is used for the following analysis.
#final model
chosen_model <- "m5"
Model <- Models$model[Models$name == chosen_model][[1]]
summary(Model)Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite
quadrature approximation with 10 quadrature points
formula: wearlikelihoodrating_home ~ wearing_position + (1 | record_id)
data: data
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 1160 -1861.86 3751.72 1483(8824) 2.94e-03 2.0e+02
Random effects:
Groups Name Variance Std.Dev.
record_id (Intercept) 3.241 1.8
Number of groups: record_id 145
Coefficients:
Estimate Std. Error z value Pr(>|z|)
wearing_positionWrist -0.7833 0.2146 -3.650 0.000263 ***
wearing_positionNecklace -0.6834 0.2138 -3.197 0.001390 **
wearing_positionSleeve collar -1.5234 0.2205 -6.908 4.92e-12 ***
wearing_positionCollar pin -1.4878 0.2170 -6.857 7.04e-12 ***
wearing_positionNeck loop -1.9235 0.2213 -8.691 < 2e-16 ***
wearing_positionHat pin -3.1911 0.2339 -13.645 < 2e-16 ***
wearing_positionGlasses -2.6719 0.2275 -11.746 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Threshold coefficients:
Estimate Std. Error z value
Extremely unlikely|Very Unlikely -4.4846 0.2543 -17.634
Very Unlikely|Unlikely -3.7356 0.2432 -15.358
Unlikely|Neutral -2.5400 0.2303 -11.028
Neutral|Likely -1.4360 0.2225 -6.454
Likely|Very likely 0.4905 0.2180 2.250
Very likely|Extremely likely 2.3480 0.2383 9.855#McFadden
McFadden(Model, Models$model[[6]])[1] 0.07320812#Significance matrix
sig_matrix <- data_significance_matrix(Model)
P4 <- Significance_matrix(sig_matrix)
ggsave(glue("output/04_unused/significance_matrix_{parameter}.pdf"),
P4, width = 4, height = 4.7)
P4
#if one wanted to print predictions for different subject percentiles
stDev <- as.numeric(attr(VarCorr(Model)$record_id, "stddev"))^2
# 5th percentile:
prob_table(Model$beta * qnorm(0.05) * stDev,
"Predicted probabilities for the **5th percentile** for **{parameter}**")| Wearing Position | Extremely unlikely | Very Unlikely | Unlikely | Neutral | Likely | Very likely | Extremely likely |
|---|---|---|---|---|---|---|---|
| Chest Pin | 1.12% | 1.21% | 4.98% | 11.91% | 42.81% | 29.26% | 8.72% |
| Wrist | 0.02% | 0.02% | 0.08% | 0.24% | 2.08% | 11.40% | 86.16% |
| Necklace | 0.03% | 0.03% | 0.14% | 0.41% | 3.48% | 17.39% | 78.51% |
| Sleeve collar | 0.00% | 0.00% | 0.00% | 0.00% | 0.04% | 0.26% | 99.69% |
| Collar pin | 0.00% | 0.00% | 0.00% | 0.01% | 0.05% | 0.32% | 99.63% |
| Neck loop | 0.00% | 0.00% | 0.00% | 0.00% | 0.00% | 0.03% | 99.96% |
| Hat pin | 0.00% | 0.00% | 0.00% | 0.00% | 0.00% | 0.00% | 100.00% |
| Glasses | 0.00% | 0.00% | 0.00% | 0.00% | 0.00% | 0.00% | 100.00% |
#average:
prob_table(Model$beta,
"Predicted probabilities for the **average person** for **{parameter}**")| Wearing Position | Extremely unlikely | Very Unlikely | Unlikely | Neutral | Likely | Very likely | Extremely likely |
|---|---|---|---|---|---|---|---|
| Chest Pin | 1.12% | 1.21% | 4.98% | 11.91% | 42.81% | 29.26% | 8.72% |
| Wrist | 2.41% | 2.55% | 9.76% | 19.52% | 43.90% | 17.68% | 4.18% |
| Necklace | 2.19% | 2.33% | 9.00% | 18.52% | 44.36% | 19.01% | 4.60% |
| Sleeve collar | 4.92% | 4.95% | 16.70% | 25.62% | 36.04% | 9.73% | 2.04% |
| Collar pin | 4.76% | 4.80% | 16.33% | 25.42% | 36.55% | 10.04% | 2.11% |
| Neck loop | 7.17% | 6.87% | 21.02% | 26.89% | 29.84% | 6.83% | 1.38% |
| Hat pin | 21.53% | 15.19% | 29.01% | 19.53% | 12.28% | 2.06% | 0.39% |
| Glasses | 14.03% | 11.63% | 27.63% | 24.19% | 18.45% | 3.40% | 0.66% |
# 95th percentile:
prob_table(Model$beta * qnorm(0.95) * stDev,
"Predicted probabilities for the **95th percentile** for **{parameter}**")| Wearing Position | Extremely unlikely | Very Unlikely | Unlikely | Neutral | Likely | Very likely | Extremely likely |
|---|---|---|---|---|---|---|---|
| Chest Pin | 1.12% | 1.21% | 4.98% | 11.91% | 42.81% | 29.26% | 8.72% |
| Wrist | 42.35% | 18.49% | 22.86% | 10.23% | 5.13% | 0.78% | 0.15% |
| Necklace | 30.13% | 17.57% | 27.39% | 15.00% | 8.33% | 1.33% | 0.25% |
| Sleeve collar | 97.44% | 1.34% | 0.85% | 0.25% | 0.11% | 0.02% | 0.00% |
| Collar pin | 96.92% | 1.60% | 1.03% | 0.30% | 0.13% | 0.02% | 0.00% |
| Neck loop | 99.69% | 0.16% | 0.10% | 0.03% | 0.01% | 0.00% | 0.00% |
| Hat pin | 100.00% | 0.00% | 0.00% | 0.00% | 0.00% | 0.00% | 0.00% |
| Glasses | 99.99% | 0.00% | 0.00% | 0.00% | 0.00% | 0.00% | 0.00% |
#random effect distribution with random effects
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:145, labels=ord.re)
axis(2)
for(i in 1:145) segments(i, ci[i,1], i, ci[i,2])
abline(h=0, lty=2)
#without confidence intervals
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:145, labels=ord.re)
axis(2)
abline(h=0, lty=2)
parameter <- "wearlikelihoodrating_social"
parameter_symbol <- sym(parameter)
#model generation
Models <- Infe(!!parameter_symbol, data_long)
#model comparisons
Model_Comparisons(Models) %>%
gt(caption = md(glue("Model comparisons for **{parameter}**"))) %>%
fmt_number(decimals = 3)| desc | p.value |
|---|---|
| interaction of sex and wearing position | 0.032 |
| inclusion of sex | 0.824 |
| inclusion of sample location | 0.129 |
| inclusion of wearing position | <0.001 |
The data supports Wearing Position, and its interaction with Sex. Sample is not supported by the data. The resulting model m4 is used for the following analysis.
#final model
chosen_model <- "m4"
Model <- Models$model[Models$name == chosen_model][[1]]
summary(Model)Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite
quadrature approximation with 10 quadrature points
formula: wearlikelihoodrating_social ~ wearing_position * demographics_sex +
(1 | record_id)
data: data
subset: demographics_sex != "Other"
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 1152 -1810.43 3664.86 2853(19418) 1.96e-03 6.7e+02
Random effects:
Groups Name Variance Std.Dev.
record_id (Intercept) 3.094 1.759
Number of groups: record_id 144
Coefficients:
Estimate Std. Error z value
wearing_positionWrist -0.86456 0.31040 -2.785
wearing_positionNecklace -1.39967 0.30101 -4.650
wearing_positionSleeve collar -2.67891 0.30999 -8.642
wearing_positionCollar pin -2.25883 0.31001 -7.286
wearing_positionNeck loop -2.84337 0.31532 -9.017
wearing_positionHat pin -3.36561 0.31872 -10.560
wearing_positionGlasses -3.63839 0.33821 -10.758
demographics_sexMale -0.61176 0.41761 -1.465
wearing_positionWrist:demographics_sexMale 0.66775 0.43052 1.551
wearing_positionNecklace:demographics_sexMale 0.70436 0.42624 1.652
wearing_positionSleeve collar:demographics_sexMale 1.20096 0.43381 2.768
wearing_positionCollar pin:demographics_sexMale 0.02352 0.42906 0.055
wearing_positionNeck loop:demographics_sexMale 0.98738 0.43120 2.290
wearing_positionHat pin:demographics_sexMale 1.10079 0.43432 2.535
wearing_positionGlasses:demographics_sexMale 1.15958 0.45243 2.563
Pr(>|z|)
wearing_positionWrist 0.00535 **
wearing_positionNecklace 3.32e-06 ***
wearing_positionSleeve collar < 2e-16 ***
wearing_positionCollar pin 3.19e-13 ***
wearing_positionNeck loop < 2e-16 ***
wearing_positionHat pin < 2e-16 ***
wearing_positionGlasses < 2e-16 ***
demographics_sexMale 0.14295
wearing_positionWrist:demographics_sexMale 0.12090
wearing_positionNecklace:demographics_sexMale 0.09844 .
wearing_positionSleeve collar:demographics_sexMale 0.00563 **
wearing_positionCollar pin:demographics_sexMale 0.95629
wearing_positionNeck loop:demographics_sexMale 0.02203 *
wearing_positionHat pin:demographics_sexMale 0.01126 *
wearing_positionGlasses:demographics_sexMale 0.01038 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Threshold coefficients:
Estimate Std. Error z value
Extremely unlikely|Very Unlikely -4.0338 0.3200 -12.607
Very Unlikely|Unlikely -2.9256 0.3105 -9.423
Unlikely|Neutral -1.7688 0.3031 -5.836
Neutral|Likely -0.6461 0.2987 -2.163
Likely|Very likely 1.2000 0.3017 3.978
Very likely|Extremely likely 3.2553 0.3490 9.329#McFadden
McFadden(Model, Models$model[[6]])[1] 0.08991832#Significance matrix
sig_matrix <- data_significance_matrix(Model, subset = "true")
sig_matrix <-
sig_matrix %>% mutate(term = str_remove(term, "demographics_sex"))
P5 <-
Significance_matrix(
sig_matrix %>% filter(!str_detect(term, "Male")), sex = TRUE
) +
geom_tile(
data = sig_matrix %>% filter(str_detect(term,"Male") & different) %>%
mutate(term = str_remove(term, ":Male")),
fill = "red")
ggsave(glue("output/04_unused/significance_matrix_{parameter}.pdf"),
P5, width = 4, height = 4.7)
P5 
#random effect distribution with random effects
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:144, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:144, labels=ord.re)
axis(2)
for(i in 1:144) segments(i, ci[i,1], i, ci[i,2])
abline(h=0, lty=2)
#without confidence intervals
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:144, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:144, labels=ord.re)
axis(2)
abline(h=0, lty=2)
parameter <- "wearlikelihoodrating_exercise"
parameter_symbol <- sym(parameter)
#model generation
Models <- Infe(!!parameter_symbol, data_long)
#model comparisons
Model_Comparisons(Models) %>%
gt(caption = md(glue("Model comparisons for **{parameter}**"))) %>%
fmt_number(decimals = 3)| desc | p.value |
|---|---|
| interaction of sex and wearing position | 0.051 |
| inclusion of sex | 0.61 |
| inclusion of sample location | 0.506 |
| inclusion of wearing position | <0.001 |
The only parameter which is supported by the data is Wearing Position. Sex, its interaction with Wearing Position, and Sample are not supported by the data. The resulting model m5 is used for the following analysis.
#final model
chosen_model <- "m5"
Model <- Models$model[Models$name == chosen_model][[1]]
summary(Model)Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite
quadrature approximation with 10 quadrature points
formula: wearlikelihoodrating_exercise ~ wearing_position + (1 | record_id)
data: data
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 1160 -1915.38 3858.76 1662(8332) 1.14e-03 2.2e+02
Random effects:
Groups Name Variance Std.Dev.
record_id (Intercept) 1.936 1.392
Number of groups: record_id 145
Coefficients:
Estimate Std. Error z value Pr(>|z|)
wearing_positionWrist 0.08459 0.21428 0.395 0.693
wearing_positionNecklace -1.97435 0.21822 -9.048 < 2e-16 ***
wearing_positionSleeve collar -0.12314 0.21650 -0.569 0.570
wearing_positionCollar pin -1.22176 0.21253 -5.749 8.99e-09 ***
wearing_positionNeck loop -2.01988 0.21945 -9.204 < 2e-16 ***
wearing_positionHat pin -2.52201 0.22273 -11.323 < 2e-16 ***
wearing_positionGlasses -3.12963 0.23294 -13.435 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Threshold coefficients:
Estimate Std. Error z value
Extremely unlikely|Very Unlikely -3.6177 0.2210 -16.369
Very Unlikely|Unlikely -2.7594 0.2108 -13.088
Unlikely|Neutral -1.6150 0.2003 -8.062
Neutral|Likely -0.7049 0.1947 -3.620
Likely|Very likely 0.9452 0.1962 4.817
Very likely|Extremely likely 2.6629 0.2271 11.726#McFadden
McFadden(Model, Models$model[[6]])[1] 0.09162989#Significance matrix
sig_matrix <- data_significance_matrix(Model)
P6 <- Significance_matrix(sig_matrix)
ggsave(glue("output/04_unused/significance_matrix_{parameter}.pdf"),
P6, width = 4, height = 4.7)
P6
#if one wanted to print predictions for different subject percentiles
stDev <- as.numeric(attr(VarCorr(Model)$record_id, "stddev"))^2
# 5th percentile:
prob_table(Model$beta * qnorm(0.05) * stDev,
"Predicted probabilities for the **5th percentile** for **{parameter}**")| Wearing Position | Extremely unlikely | Very Unlikely | Unlikely | Neutral | Likely | Very likely | Extremely likely |
|---|---|---|---|---|---|---|---|
| Chest Pin | 2.61% | 3.34% | 10.63% | 16.48% | 38.94% | 21.46% | 6.52% |
| Wrist | 3.40% | 4.26% | 13.00% | 18.62% | 37.83% | 17.83% | 5.06% |
| Necklace | 0.00% | 0.01% | 0.03% | 0.05% | 0.38% | 2.12% | 97.41% |
| Sleeve collar | 1.78% | 2.32% | 7.74% | 13.18% | 38.45% | 27.16% | 9.36% |
| Collar pin | 0.05% | 0.07% | 0.28% | 0.59% | 3.99% | 17.65% | 77.36% |
| Neck loop | 0.00% | 0.01% | 0.02% | 0.05% | 0.33% | 1.84% | 97.75% |
| Hat pin | 0.00% | 0.00% | 0.00% | 0.01% | 0.07% | 0.38% | 99.54% |
| Glasses | 0.00% | 0.00% | 0.00% | 0.00% | 0.01% | 0.06% | 99.93% |
#average:
prob_table(Model$beta,
"Predicted probabilities for the **average person** for **{parameter}**")| Wearing Position | Extremely unlikely | Very Unlikely | Unlikely | Neutral | Likely | Very likely | Extremely likely |
|---|---|---|---|---|---|---|---|
| Chest Pin | 2.61% | 3.34% | 10.63% | 16.48% | 38.94% | 21.46% | 6.52% |
| Wrist | 2.41% | 3.09% | 9.95% | 15.78% | 39.05% | 22.67% | 7.06% |
| Necklace | 16.20% | 15.12% | 27.56% | 19.18% | 16.82% | 4.16% | 0.96% |
| Sleeve collar | 2.95% | 3.74% | 11.68% | 17.49% | 38.58% | 19.76% | 5.81% |
| Collar pin | 8.35% | 9.34% | 22.60% | 22.35% | 27.08% | 8.26% | 2.01% |
| Neck loop | 16.83% | 15.48% | 27.67% | 18.85% | 16.26% | 3.99% | 0.92% |
| Hat pin | 25.05% | 19.04% | 27.14% | 14.78% | 10.95% | 2.47% | 0.56% |
| Glasses | 38.04% | 21.12% | 22.82% | 9.90% | 6.46% | 1.37% | 0.30% |
# 95th percentile:
prob_table(Model$beta * qnorm(0.95) * stDev,
"Predicted probabilities for the **95th percentile** for **{parameter}**")| Wearing Position | Extremely unlikely | Very Unlikely | Unlikely | Neutral | Likely | Very likely | Extremely likely |
|---|---|---|---|---|---|---|---|
| Chest Pin | 2.61% | 3.34% | 10.63% | 16.48% | 38.94% | 21.46% | 6.52% |
| Wrist | 2.01% | 2.60% | 8.57% | 14.21% | 38.88% | 25.35% | 8.37% |
| Necklace | 93.53% | 3.62% | 1.92% | 0.55% | 0.30% | 0.06% | 0.01% |
| Sleeve collar | 3.82% | 4.75% | 14.17% | 19.50% | 36.96% | 16.29% | 4.50% |
| Collar pin | 56.80% | 18.82% | 15.07% | 5.34% | 3.18% | 0.64% | 0.14% |
| Neck loop | 94.35% | 3.17% | 1.67% | 0.48% | 0.26% | 0.05% | 0.01% |
| Hat pin | 98.80% | 0.68% | 0.35% | 0.10% | 0.05% | 0.01% | 0.00% |
| Glasses | 99.83% | 0.10% | 0.05% | 0.01% | 0.01% | 0.00% | 0.00% |
#random effect distribution with random effects
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:145, labels=ord.re)
axis(2)
for(i in 1:145) segments(i, ci[i,1], i, ci[i,2])
abline(h=0, lty=2)
#without confidence intervals
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:145, labels=ord.re)
axis(2)
abline(h=0, lty=2)
parameter <- "wearduration"
parameter_symbol <- sym(parameter)
#model generation
Models <- Infe(!!parameter_symbol, data_long)
#model comparisons
Model_Comparisons(Models) %>%
gt(caption = md(glue("Model comparisons for **{parameter}**"))) %>%
fmt_number(decimals = 3)| desc | p.value |
|---|---|
| interaction of sex and wearing position | 0.005 |
| inclusion of sex | 0.852 |
| inclusion of sample location | 0.852 |
| inclusion of wearing position | <0.001 |
The data supports Wearing Position, and its interaction with Sex. Sample is not supported by the data. The resulting model m4 is used for the following analysis.
#final model
chosen_model <- "m4"
Model <- Models$model[Models$name == chosen_model][[1]]
summary(Model)Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite
quadrature approximation with 10 quadrature points
formula: wearduration ~ wearing_position * demographics_sex + (1 | record_id)
data: data
subset: demographics_sex != "Other"
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 1152 -1904.65 3853.30 3096(16129) 1.60e-03 7.1e+02
Random effects:
Groups Name Variance Std.Dev.
record_id (Intercept) 2.521 1.588
Number of groups: record_id 144
Coefficients:
Estimate Std. Error z value
wearing_positionWrist -0.7062 0.3044 -2.320
wearing_positionNecklace -1.4241 0.2999 -4.749
wearing_positionSleeve collar -2.3526 0.3103 -7.582
wearing_positionCollar pin -2.0286 0.3094 -6.557
wearing_positionNeck loop -3.2638 0.3185 -10.247
wearing_positionHat pin -4.0880 0.3266 -12.518
wearing_positionGlasses -3.8876 0.3300 -11.779
demographics_sexMale -0.8566 0.4047 -2.117
wearing_positionWrist:demographics_sexMale 0.5400 0.4336 1.245
wearing_positionNecklace:demographics_sexMale 0.3419 0.4268 0.801
wearing_positionSleeve collar:demographics_sexMale 1.1391 0.4337 2.627
wearing_positionCollar pin:demographics_sexMale 0.3680 0.4317 0.853
wearing_positionNeck loop:demographics_sexMale 1.3336 0.4350 3.066
wearing_positionHat pin:demographics_sexMale 1.5122 0.4376 3.456
wearing_positionGlasses:demographics_sexMale 1.1499 0.4450 2.584
Pr(>|z|)
wearing_positionWrist 0.020348 *
wearing_positionNecklace 2.04e-06 ***
wearing_positionSleeve collar 3.39e-14 ***
wearing_positionCollar pin 5.48e-11 ***
wearing_positionNeck loop < 2e-16 ***
wearing_positionHat pin < 2e-16 ***
wearing_positionGlasses < 2e-16 ***
demographics_sexMale 0.034286 *
wearing_positionWrist:demographics_sexMale 0.212968
wearing_positionNecklace:demographics_sexMale 0.423084
wearing_positionSleeve collar:demographics_sexMale 0.008626 **
wearing_positionCollar pin:demographics_sexMale 0.393928
wearing_positionNeck loop:demographics_sexMale 0.002171 **
wearing_positionHat pin:demographics_sexMale 0.000549 ***
wearing_positionGlasses:demographics_sexMale 0.009758 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Threshold coefficients:
Estimate Std. Error z value
Not at all|Up to 1 hour -4.6898 0.3164 -14.820
Up to 1 hour|1-3 hours -3.6121 0.3061 -11.802
1-3 hours|3-6 hours -2.2042 0.2950 -7.471
3-6 hours|6-8 hours -1.0963 0.2886 -3.799
6-8 hours|8-12 hours -0.1181 0.2855 -0.414
8-12 hours|Throughout the day 0.6863 0.2871 2.390#McFadden
McFadden(Model, Models$model[[6]])[1] 0.1005441#Significance matrix
sig_matrix <- data_significance_matrix(Model, subset = "true")
sig_matrix <-
sig_matrix %>% mutate(term = str_remove(term, "demographics_sex"))
P7 <-
Significance_matrix(
sig_matrix %>% filter(!str_detect(term, "Male")), sex = TRUE
) +
geom_tile(
data = sig_matrix %>% filter(str_detect(term,"Male") & different) %>%
mutate(term = str_remove(term, ":Male")),
fill = "red")
ggsave(glue("output/04_unused/significance_matrix_{parameter}.pdf"),
P7, width = 4, height = 4.7)
P7
#random effect distribution with random effects
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:144, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:144, labels=ord.re)
axis(2)
for(i in 1:144) segments(i, ci[i,1], i, ci[i,2])
abline(h=0, lty=2)
#without confidence intervals
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:144, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:144, labels=ord.re)
axis(2)
abline(h=0, lty=2)
parameter <- "restrict"
parameter_symbol <- sym(parameter)
#model generation
Models <- Infe(!!parameter_symbol, data_long)
#model comparisons
Model_Comparisons(Models) %>%
gt(caption = md(glue("Model comparisons for **{parameter}**"))) %>%
fmt_number(decimals = 3)| desc | p.value |
|---|---|
| interaction of sex and wearing position | 0.3 |
| inclusion of sex | 0.918 |
| inclusion of sample location | 0.918 |
| inclusion of wearing position | <0.001 |
The only parameter which is supported by the data is Wearing Position. Sex, its interaction with Wearing Position, and Sample are not supported by the data. The resulting model m5 is used for the following analysis.
#final model
chosen_model <- "m5"
Model <- Models$model[Models$name == chosen_model][[1]]
summary(Model)Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite
quadrature approximation with 10 quadrature points
formula: restrict ~ wearing_position + (1 | record_id)
data: data
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 1160 -1764.02 3556.03 1587(4764) 9.38e-04 3.4e+02
Random effects:
Groups Name Variance Std.Dev.
record_id (Intercept) 1.201 1.096
Number of groups: record_id 145
Coefficients:
Estimate Std. Error z value Pr(>|z|)
wearing_positionWrist -1.3797 0.2457 -5.614 1.97e-08 ***
wearing_positionNecklace -1.9653 0.2436 -8.069 7.11e-16 ***
wearing_positionSleeve collar -1.4102 0.2523 -5.590 2.27e-08 ***
wearing_positionCollar pin -1.9734 0.2489 -7.927 2.25e-15 ***
wearing_positionNeck loop -2.3979 0.2478 -9.678 < 2e-16 ***
wearing_positionHat pin -2.6988 0.2491 -10.836 < 2e-16 ***
wearing_positionGlasses -3.1695 0.2561 -12.374 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Threshold coefficients:
Estimate Std. Error z value
Extremely|Significantly -6.3272 0.3051 -20.735
Significantly|Quite a bit -5.2148 0.2633 -19.807
Quite a bit|Moderately -4.2970 0.2450 -17.540
Moderately|Somewhat -3.4811 0.2335 -14.908
Somewhat|Slightly -2.4176 0.2222 -10.882
Slightly|Not at all -0.9715 0.2107 -4.612#McFadden
McFadden(Model, Models$model[[6]])[1] 0.06010535#Significance matrix
sig_matrix <- data_significance_matrix(Model)
P8 <- Significance_matrix(sig_matrix)
ggsave(glue("output/04_unused/significance_matrix_{parameter}.pdf"),
P8, width = 4, height = 4.7)
P8
#if one wanted to print predictions for different subject percentiles
stDev <- as.numeric(attr(VarCorr(Model)$record_id, "stddev"))^2
# 5th percentile:
prob_table(Model$beta * qnorm(0.05) * stDev,
"Predicted probabilities for the **5th percentile** for **{parameter}**")| Wearing Position | Extremely | Significantly | Quite a bit | Moderately | Somewhat | Slightly | Not at all |
|---|---|---|---|---|---|---|---|
| Chest Pin | 0.18% | 0.36% | 0.80% | 1.64% | 5.20% | 19.27% | 72.54% |
| Wrist | 0.01% | 0.02% | 0.05% | 0.11% | 0.38% | 1.84% | 97.58% |
| Necklace | 0.00% | 0.01% | 0.02% | 0.04% | 0.12% | 0.59% | 99.23% |
| Sleeve collar | 0.01% | 0.02% | 0.05% | 0.11% | 0.36% | 1.74% | 97.72% |
| Collar pin | 0.00% | 0.01% | 0.02% | 0.03% | 0.12% | 0.58% | 99.24% |
| Neck loop | 0.00% | 0.00% | 0.01% | 0.02% | 0.05% | 0.25% | 99.67% |
| Hat pin | 0.00% | 0.00% | 0.00% | 0.01% | 0.03% | 0.14% | 99.82% |
| Glasses | 0.00% | 0.00% | 0.00% | 0.00% | 0.01% | 0.06% | 99.93% |
#average:
prob_table(Model$beta,
"Predicted probabilities for the **average person** for **{parameter}**")| Wearing Position | Extremely | Significantly | Quite a bit | Moderately | Somewhat | Slightly | Not at all |
|---|---|---|---|---|---|---|---|
| Chest Pin | 0.18% | 0.36% | 0.80% | 1.64% | 5.20% | 19.27% | 72.54% |
| Wrist | 0.71% | 1.41% | 3.02% | 5.77% | 15.26% | 33.91% | 39.94% |
| Necklace | 1.26% | 2.48% | 5.12% | 9.16% | 20.87% | 34.10% | 27.02% |
| Sleeve collar | 0.73% | 1.45% | 3.10% | 5.91% | 15.55% | 34.05% | 39.21% |
| Collar pin | 1.27% | 2.49% | 5.15% | 9.21% | 20.95% | 34.07% | 26.86% |
| Neck loop | 1.93% | 3.71% | 7.38% | 12.27% | 24.22% | 31.13% | 19.37% |
| Hat pin | 2.59% | 4.89% | 9.35% | 14.56% | 25.60% | 27.92% | 15.09% |
| Glasses | 4.08% | 7.37% | 13.01% | 17.81% | 25.69% | 22.05% | 9.99% |
# 95th percentile:
prob_table(Model$beta * qnorm(0.95) * stDev,
"Predicted probabilities for the **95th percentile** for **{parameter}**")| Wearing Position | Extremely | Significantly | Quite a bit | Moderately | Somewhat | Slightly | Not at all |
|---|---|---|---|---|---|---|---|
| Chest Pin | 0.18% | 0.36% | 0.80% | 1.64% | 5.20% | 19.27% | 72.54% |
| Wrist | 2.65% | 5.00% | 9.54% | 14.76% | 25.68% | 27.61% | 14.76% |
| Necklace | 7.98% | 12.89% | 18.90% | 20.12% | 21.33% | 13.62% | 5.16% |
| Sleeve collar | 2.81% | 5.28% | 9.97% | 15.21% | 25.82% | 26.89% | 14.02% |
| Collar pin | 8.10% | 13.04% | 19.02% | 20.12% | 21.19% | 13.45% | 5.09% |
| Neck loop | 16.93% | 21.34% | 22.55% | 17.01% | 13.22% | 6.69% | 2.26% |
| Hat pin | 26.96% | 25.93% | 20.87% | 12.64% | 8.44% | 3.89% | 1.26% |
| Glasses | 48.33% | 25.66% | 13.69% | 6.46% | 3.75% | 1.60% | 0.50% |
#random effect distribution with random effects
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:145, labels=ord.re)
axis(2)
for(i in 1:145) segments(i, ci[i,1], i, ci[i,2])
abline(h=0, lty=2)
#without confidence intervals
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:145, labels=ord.re)
axis(2)
abline(h=0, lty=2)
a <- theme(plot.margin = margin(5, 10,5,5,"pt"))
b <- theme(plot.margin = margin(5, 5,5,10,"pt"))
(P1 + a + P2 + b) / (P3 + a + P4 + b) / (P5 + a + P6 + b) / (P7 + a + P8 + b) +
plot_annotation(tag_levels = "A") & theme(plot.tag = element_text(size = 20))
ggsave("output/01_figures/figure2_summary_graph_differences.pdf", width = 10, height = 20)
ggsave("output/01_figures/figure2_summary_graph_differences.jpeg", width = 10, height = 20, dpi = 300)